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Abstract 

In the implicit large eddy simulation (ILES) paradigm, the dissipative nature of high-resolution 
shock-capturing schemes is exploited to provide an implicit model of turbulence. The ILES approach has 
been applied to different contexts, with varying degrees of success. It is the de-facto standard in many 
astrophysical simulations and in particular in studies of core-collapse supernovae (CCSN). Recent 3D 
simulations suggest that turbulence might play a crucial role in core-collapse supernova explosions, 
however the fidelity with which turbulence is simulated in these studies is unclear. Especially considering 
that the accuracy of ILES for the regime of interest in CCSN, weakly compressible and strongly 
anisotropic, has not been systematically assessed before. Anisotropy, in particular, could impact the 
dissipative properties of the flow and enhance the turbulent pressure in the radial direction, favouring the 
explosion. In this paper we assess the accuracy of ILES using numerical methods most commonly 
employed in computational astrophysics by means of a number of local simulations of driven, weakly 
compressible, anisotropic turbulence. Our simulations employ several different methods and span a wide 
range of resolutions. We report a detailed analysis of the way in which the turbulent cascade is 
influenced by the numerics. Our results suggest that anisotropy and compressibility in CCSN turbulence 
have little effect on the turbulent kinetic energy spectrum and a Kolmogorov k~ 5/3 scaling is obtained in 
the inertial range. We find that, on the one hand, the kinetic energy dissipation rate at large scales is 
correctly captured even at low resolutions, suggesting that very high “effective Reynolds number” can be 
achieved at the largest scales of the simulation. On the other hand, the dynamics at intermediate scales 
appears to be completely dominated by the so-called bottleneck effect, i.e., the pile up of kinetic energy 
close to the dissipation range due to the partial suppression of the energy cascade by numerical 
viscosity. An inertial range is not recovered until the point where high resolution ~ 512 3 , which would be 
difficult to realize in global simulations, is reached. We discuss the consequences for CCSN simulations. 
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1 Introduction 

Despite decades of studies and compelling evidence that a 
significant fraction [1] of stars with initial masses in excess 
of ~8 solar masses explode as core-collapse supernovae 
(CCSN) at the end of their evolution, the exact details of the 
explosion mechanism are still uncertain [2, 3,4, 5]. Current 
state-of-the art 3D simulations either fail to explode or have 
explosion energies that fall short of the observed energies 
by factors of a few for most of the progenitor mass range 
[6,4,5], 

The dynamics at the center of a star undergoing core col¬ 
lapse is shaped by a delicate balance between competing 

Full list of author information is available at the end of the article 


effects where all of the known forces: gravity, electromag¬ 
netism, weak and strong interactions, are important. The 
task of modeling these systems is made particularly chal¬ 
lenging by the fact that the generation of the asymptotic 
explosion energies, although enormous (~ 10 44 J), requires 
a rather subtle, percent-level imbalance between non-linear 
processes over many dynamical times. 

The flow of plasma in the core of a star going supernova 
is known to be unstable to convection [7, 8, 9, 10] and/or 
to another large scale instability known as standing accre¬ 
tion shock instability [11, 12], In any case, given the very 
large Reynolds numbers, as large as ~ 10 1 ' in the region 
of interest [13] (the so-called gain region, where neutrino 
heating dominates over neutrino cooling), it is expected 
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that the resulting flow will be fully turbulent. It has been 
suggested [14, 15] recently that turbulence and, in partic¬ 
ular, turbulent pressure could tip the balance of the forces 
in favor of explosion. In this respect, anisotropy is of key 
importance, because it results in an effective radial pres¬ 
sure support with adiabatic index 7 t ur b = 2 , much larger 
than that of thermal (radiation) pressure ( 7 ^ ~ 4/3). This 
means that turbulent kinetic energy is a much more valu¬ 
able source of radial pressure support than thermal energy 
(see Appendix A). 

All of the current numerical simulations employ the 
implicit large eddy simulation (ILES) paradigm [16, 17] 
(also known as monotone integrated LES (MILES)) of ex¬ 
ploiting the dissipative nature of high resolution shock cap¬ 
turing (HRSC) methods as an implicit turbulence model. 
However, the combination of the use of rather dissipative 
schemes and the relatively low spatial resolution that can 
be achieved in global simulations is such that the fidelity 
with which turbulence is captured is questionable [13]. 

To be useful in the context of CCSN simulations, an ILES 
should, at the very least, account for the right rate of decay 
of the kinetic energy at the largest scales while avoiding 
unphysical pile up of energy at smaller scales. Unfortu¬ 
nately, all of the current simulations seem to be strongly 
dominated by the so-called bottleneck effect [13], which 
corresponds to an inefficient energy transfer across inter¬ 
mediate scales due to the viscous suppression of non-linear 
interaction with smaller scales [18, 19, 20, 21, 22]. Current 
global simulations achieve resolutions, in the turbulent re¬ 
gion, comparable to those of 30 3 — 70 3 lattices in periodic 
domains [23, 15, 13]. At these resolutions, almost all of the 
dynamical range of the simulations can be expected to be 
directly affected by numerical viscosity [24]. The fidelity 
with which turbulence is captured in these simulations will 
then depend on the degree with which the numerical trun¬ 
cation error approximates an LES closure. 

In this respect, it has been shown by [25] and [26] that 
many HRSC methods can be too dissipative to yield a 
faithful description of turbulence at low resolutions. These 
studies, however, considered a different regime, decaying 
isotropic turbulence, while turbulence in a core-collapse su¬ 
pernova, as well as in many other astrophysical settings, is 
often strongly anisotropic [27, 14, 15] as rotational invari¬ 
ance is broken by gravity. [25] and [26] also considered 
different numerical schemes with respect to those used in 
supernova simulations. Both of these aspects can, in princi¬ 
ple, be important. First of all, strong anisotropies could po¬ 
tentially influence the turbulence dynamics at the level of 
the energy cascade and of the dissipation [28]. Secondly, 
some of the schemes used in computational astrophysics, 
such as the piecewise parabolic method (PPM) [29] as well 
as some of the MUSCL [30] schemes, have been shown, 
differently from some of the methods considered by [25] 
and [26], to be well suited for ILES [31, 32], 


The aim of this work is to fill the gap between existing 
theoretical studies and the particular applications of our in¬ 
terest. To this end we use a publicly available code, FLASH 
[33, 34, 35], which is widely used in the computational 
astrophysics community, and perform a series of simula¬ 
tions of turbulence in a regime relevant for core-collapse 
supernovae: driven at large scale, with large anisotropies 
and mildly compressible. We use five different numerical 
setups and, for each, several resolutions in the range from 
64 3 to 512 3 in a periodic domain. We study in detail the 
way in which the energy cascade across different scales is 
represented by our ILES and we discuss the use of local 
or lower dimensional diagnostics that can be used to assess 
the quality of a global simulation in a complex geometry 
where 3D spectra are not readily available. 

The rest of this paper is organized as follows. First, in 
Section 2, we discuss the exact setup of our simulations 
and the diagnostic quantities used in our analysis. Then, in 
Section 3, we discuss the basic characteristics of the flow 
realized in our simulations. In Section 4, we present a de¬ 
tailed analysis of the way in which the energy cascade is 
captured by the different schemes at different scales. In 
particular, we quantify the accuracy with which different 
methods capture the decay rate of energy from the largest 
scales and the way in which energy is distributed across 
scales. We discuss the role of anisotropies in the context 
of the 4/5—law, a fundamental exact relation for isotropic 
and incompressible turbulence relating the statistics of ve¬ 
locity fluctuations with the energy dissipation rate (see Sec¬ 
tion 2.3), in Section 5. We explore the use of the 2D, trans¬ 
verse, energy spectrum as a diagnostic for 3D simulations 
in Section 6 . Finally, we present a brief summary of our 
main findings, as well as a discussion of their implications 
for CCSN simulations in Section 7. Appendix A, contains 
some supplemental background material on the role of tur¬ 
bulence in the explosion mechanism of CCSN. 

2 Methods 

2.1 Numerical methods 

We consider a compressible fluid with a prescribed accel¬ 
eration, a, in a unit-box with periodic boundary conditions. 
The code that we employ for these simulations, FLASH, 
solves the gas-dynamics equations in conservation form. In 
particular we evolve the continuity equation 

d t p+\7-(pv) = 0 ( 1 ) 

and the momentum equation 

(9 t (pv) + V • (pv® v+pl) = pa. (2) 

These equations are closed with a simple isentropic equa¬ 
tion of state. 


( 3 ) 
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that can be considered as a rough description of a gas dom¬ 
inated by radiation pressure. Since the equation of state en¬ 
sures an adiabatic evolution we do not need to solve the en¬ 
ergy equation as equations (1), (2) and (3) suffice to fully 
describe the flow. 

Equations (1) & (2) are solved using the directionally- 
unsplit hydrodynamics solver of the open-source FLASH 
simulation framework. FLASH implements the corner 
transport upwind method [36] for fully directionally- 
unsplit evolution of the Euler equations [37, 38]. FLASH 
includes several options for the order of spatial reconstruc¬ 
tion [35], including 2 nd -order TVD [30], 3 rd -order PPM 
[29], and 5 th -order WENOZ [39]. Fluxes are computed at 
2 nd order accuracy using one of a number of approximate 
Riemann solvers included in FLASH, such as HLLE [40] 
and HLLC [41]. Second-order accuracy in time is achieved 
via a characteristic tracing evolution of the Riemann solver 
input states to the time step midpoint [29]. We remark that, 
in accordance with the ILES, paradigm, we do not include 
any additional sub-grid scale model, but relied on the im¬ 
plicit turbulent closure built in the numerical schemes we 
use for the integration of the hydrodynamics equation. 

All of our simulations start with the fluid at rest p = 1, 
v = 0. Turbulence is driven using the stirring module of 
FLASH. This module uses the Ornstein-Uhlenbeck process 
[42] to generate stirring modes in Fourier space. This yields 
an acceleration field which smoothly decorrelates [43] over 
a timescale T s . The FLASH implementation permits the 
use of any arbitrary combination of solenoidal and com¬ 
pressive modes [44]. For our runs, we set T s = 0.5, we 
use only solenoidal forcing and we restrict the accelerat¬ 
ing field to be nonzero only in the first four Fourier modes. 
This forcing is designed to mimic the influence of some 
larger scale weakly compressible flow and, for this reason, 
it does not include any compressible component. This is 
a reasonable approximation for low Mach number convec¬ 
tion which is well described by the anelastic approxima¬ 
tion, e.g., [45], In the CCSN context, simulations show that 
the turbulence is highly anisotropic, being roughly twice 
as strong in the radial direction as either tangential direc¬ 
tion [46, 14, 47, 15] since it is driven by buoyancy due to 
a negative radial entropy gradient. In order to emulate this 
behavior, the accelerating field in the direction (which 
is going to play the role of the radial direction) is scaled 
by a constant factor (before the solenoidal projection of the 
acceleration field) such that R xx ~ 2 R yy ~ 2 R zz , where 

Rij = (pvi Vj) , (4) 

is the Reynolds stress tensor (to simplify the notation we 
considered a frame in which (pv) = 0) and (•) denotes 
an ensemble average. Finally, the overall strength of the 
stirring is tuned to achieve a RMS Mach number of ~ 0.35, 
which is typically observed in realistic CCSN simulations 
[48, 49], 


2.2 Energy transfer equations 

In order to study the cascade of the specific kinetic energy 
(which we will refer to simply as “kinetic energy” or “en¬ 
ergy” in the following), |v| 2 /2, we will consider an energy 
budget equation across different scales, analogous to the 
one commonly employed in the study of incompressible, 
isotropic turbulence, e.g., [50], In particular, we consider 
the momentum equation (2) in non-conservation form, 

d t v + (v • V) v = -V Vp + a, (5) 

where V = 1/p is the specific volume of the gas. 

We can use equation (5) to derive an evolution equation 


for the Fourier transform of the velocity 

v(k) = f e- 2,rikx v(x)d 3 x. (6) 

J R 3 

Transforming both sides of equation (5) we obtain 

<9 t v + v * 27rik <8> v = —V * 27rikp + a, (7) 

where * denotes the convolution operator, i.e., 

[f*g }{ k ) = / /(q)5(k- q)d 3 q. (8) 

Jr 3 

If we multiply both sides of equation (7) by v* and take the 
real part, we obtain an equation for the 3D energy spectrum 

d t E( k) = T(k) + C(k) + e(k), (9) 

where 

E(k) = * v • v* , (10) 

T(k) = —27r 3? (v * ik (g) v) • v* , (11) 

C'(k) = —27 t (V * ikp) • v* ; (12) 

e(k) = a • v* . (13) 


Here E is the energy spectrum (the velocity power spec¬ 
tral density (PSD)) and T is the same transfer term as in 
the classical incompressible equations and e is the energy 
injection rate. The C term vanishes in the incompressible 
limit and represents the interaction between kinetic and 
acoustic modes. In practice, in our models, C is found to 
be at least one order of magnitude smaller than T at all 
scales and it is thus negligible. In any case, we retain C in 
the analysis below. 

For each of the spectral quantities, S, being E , T, C or e, 
we define the integrated spectrum, S(k), as 

S(k) = f S(k)<$(|k| - Jfe) dk, 

Jr 3 


( 14 ) 
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S(-) being the Dirac delta function. 

Integrating equation (9), we obtain the following one¬ 
dimensional energy balance equation 

d t E{k) = T{k) + C{k) + e(jfe). (15) 

This can also be written in terms of the energy flux across 
scales, 

n(fc) = - f (16) 

Jo 

as 


2.3 Structure functions 

The energy spectrum and its sources/fluxes give a compre¬ 
hensive picture of the energy cascade and can be used to 
assess the level of convergence of the simulation. Unfor¬ 
tunately, 3D energy spectra and fluxes are not easily ac¬ 
cessible in calculations in complex domains and/or with 
inhomogeneous turbulence. In these cases, local quantities 
in the physical domain are more easily extracted and ana¬ 
lyzed. Hence, one of the goals of this work is to validate the 
use of indirect measures of convergence of ILES. Among 
these quantities, the structure functions of the velocity ap¬ 
pear to be natural candidates for study. 

We define the velocity increments 


d t E{k) + d fc n(fc) = C(k) + e(k). 


(17) 


Sv(x, r) = [v(x -f r) — v(x)] 


r 

r 


( 22 ) 


Notice that we did not assume isotropy in any of the above. 

Equation (15) is derived in the inviscid limit. In practice, 
our evolution method introduces dissipation in the form of 
“numerical viscosity”. This can be quantified in terms of 
the residual 


R(k) = d t E{k) - T{k) - C{k) - e(k ). (18) 


This can be used to define a wave number dependent nu¬ 
merical viscosity: 


v(k) 


1 R(k) 

2 k 2 E(k) ' 


(19) 


We remark that v does not, in general, correspond to a clas¬ 
sical shear or bulk viscosity, but can nevertheless be inter¬ 
preted as a relative measure of the dissipation acting at dif¬ 
ferent wave numbers (see, e.g., [51, 52, 53] for alternative 
approaches). 

In practice, since we will be working in the station¬ 
ary case, after having taken the appropriate time averages, 
R(k) reduces to 


and study the quantities 

S p {r) = {Sv p )j= 0 , (23) 

where, (•}j= o denotes an ensemble average as well as a 
mean over all of the angles between v and r (in other words 
we are looking at the j = 0 component of the SO (3) de¬ 
composition of the structure functions [54]). In the case of 
homogeneous turbulence S p does not depend on x and is 
thus a function of only the separation r. 

The most important relation involving the structure func¬ 
tions is the so-called 4/5—law, which relates the third order 
structure function, £>3 (r), with the mean energy dissipation 
rate. 


(e) = / e (k) d k , 


(24) 


and states that, for incompressible, homogeneous and 
isotropic turbulence [50]: 



(25) 


R(k) = -T(k) - C(k) - e(k). (20) 

Finally, since we are working in a periodic domain, which 
we take of size L x = L y = L z = 1, all of the spectra 
are quantized and non-trivial only for k x . k y and k z inte¬ 
gers. Furthermore, all of the integrals in wave number space 
reduce to summations. Integrals over spherical shells are 
transformed to weighted sums following [43]: 

Anrt - 2 . . 

E ik) = — ^( k )’ < 21 ) 

k fc-l/2<|k|<fc+l/2 

where is the number of discrete wave-numbers in the 
shell k- 1/2 < |fc| < k + 1 / 2 . 


Equation (25) can be derived from the Navier-Stokes equa¬ 
tion for fully-developed, incompressible, homogeneous 
and isotropic turbulence and it is one of the few exact re¬ 
lations in the theory of turbulence [50], In the anisotropic 
or compressible case, however, equation (25) is not strictly 
valid and could be violated in the data. As we show in Sec¬ 
tion 5, we find equation (25) to be very well satisfied by our 
data, suggesting that the 3 rd order structure function can be 
a very useful diagnostic in global simulations. 

2.4 Transverse energy spectrum 
Another alternative to the analysis of 3D spectra, which 
has been adopted by several authors in the core-collapse 
supernova context [55, 23, 47, 13], is the use of 2D spec¬ 
tra computed using a spherical harmonics expansion of the 
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velocity field tangential to one or more spherical shells in 
the simulation. Analogously, we emulate this by looking at 
quantities in the y — z plane and we define the 2D spectra 

E±(k_ l) = ^ J vj. -vl + kl - k±') dk y dk z , 

(26) 

where vj_ is the projection of the velocity perpendicular 
to the £—direction and we introduced the partial Fourier 
transform of vj_ : 

]_ r L *l 2 

v±(k v ,k z )= lim — dx 

L ^+oo L x J —L x /2 ^ 7 ) 

f e- 2ni ( k * v+k * z K ± (x,y,z)dydz. 

J R 2 

In the limit of infinite Reynolds number / resolution, the 2D 
spectrum is expected to have the same asymptotic behavior 
as the 3D spectrum, however it is a-priori unclear if E± is 
a good proxy for E at finite resolution. For this reason we 
find it useful to investigate this here. 

As was the case for the 3D spectra, also here the spectrum 
is non-trivial only for integer k y and k z , when periodicity is 
taken into account. The integral in equation (26) is treated 
analogously to the integral in the equation (14) for the 3D 
case, while the average in the x- direction in equation (27) 
is converted to an average over the a;—extent of the simula¬ 
tion box. 

3 Basic flow properties 

We employ the finite-volume HRSC (Godunov) approach 
in which physical states are reconstructed at inter-cell 
boundaries and local Riemann problems are solved to com¬ 
pute the physical inter cell fluxes. In particular, we per¬ 
form five groups of simulations using different numerical 
methods. Each group is labeled using the name of the re¬ 
construction algorithm and of the Riemann solver. For in¬ 
stance TVD.HLLE, denotes a group of simulations done us¬ 
ing TVD reconstruction and HLLE Riemann solver. Single 
simulations are labeled using their resolution so that, for 
instance, TVD_HLLE_N12 8, denotes the TVD.HLLE run 
done using a 128 3 grid. For all of the runs the timestep 
is chosen to have a CFL, i.e., cAt/Ax, of 0.4, c being the 
maximum characteristics speed, with the exception of the 
PPM_HLLC_CFL0.8 runs where we set the CFL to 0.8. 
For the TVD runs we use the monotonized central (MC) 
slope limiter [30]. The runs with PPM use the original flat¬ 
tening and artificial viscosity prescriptions from [29]. The 
artificial viscosity coefficient is 0.1. We remark that the use 
of the artificial viscosity for PPM is not really necessary in 
this regime [56], however our goal is not to perform a study 
of the turbulent dynamics, but to assess how each numeri¬ 
cal method performs when used under the same condition 


as in a real CCSN simulation where strong shocks need to 
be handled in some parts of the domain. 

For each group of simulations we run four resolutions: 
64 3 ,128 3 , 256 3 and 512 3 . The RMS velocity in all of the 
runs is r) lms ~ 0.4, giving an eddy turnover time r = 
l/t 4 ms — 2.5. All of the simulations are run until time 
t = 100 (~ 40 eddy turnover times). The time evolu¬ 
tion of a few relevant diagnostics is shown in Figure 1 for 
our fiducial group of runs (PPM_HLLC) at different resolu¬ 
tions. We can see how the flow is accelerated from rest and 
quickly reaches a steady, fully turbulent, state. In all cases, 
steady state is reached after t > 3 (~ 1 turnover time) 
and the diagnostics are insensitive to the resolution. The 
results for the other runs (not shown) are very similar to the 
ones of PPM_HLLC as they all achieve very similar RMS 
Mach numbers and Reynolds stresses. All of the analysis 
shown in the rest of the paper are performed using 380 3D 
snapshots (evenly spaced in time) of the data in the interval 
5 < t < 100. 

A first, qualitative, comparison between the different 
methods can be done by looking at their visualizations. 
In particular, in Figure 2, we show a visualization of the 
magnitude of the vorticity in the x—z plane for four of the 
five schemes (excluding PPM_HLLC_CFL0.8) at the high¬ 
est resolution (512 3 ). The data is taken at the final time 
(t = 100). As it can be seen from the figure, all of the simu¬ 
lations show the presence of thin, elongated, regions of high 
vorticity, as typically seen in direct numerical simulations 
(DNS) of homogeneous turbulent flows [57, 58], How¬ 
ever, the width and the intensity of the vorticity at these 
smaller scales depend crucially on the numerical scheme. 
Methods with small intrinsic numerical viscosity, such as 
PPM_HLLC and WENOZ.HLLC, present smaller structures 
and more intermittent vorticity fields with respect to more 
dissipative methods, such as PPM.HLLE and TVD.HLLE. 

4 The energy cascade 

In this section we focus our analysis on the accuracy with 
which the energy cascade is captured by our ILES runs. 
First, we focus on the largest scales of the simulation with 
the goal of quantifying the accuracy in the decay rate of 
the energy as a function of the resolution for the different 
methods. Next, we will look at the energy distribution at 
smaller scales where, in resolved simulations, the inertial 
range starts. Finally, we will look at the dynamics in the 
dissipation region and summarize. 

4.1 Energy decay rate 

In the limit of very large Reynolds number it is assumed, 
in standard turbulence phenomenology [50], that there ex¬ 
ists a range of wave numbers (the inertial range) where en¬ 
ergy injection and dissipation can be neglected in equation 
(17). In this range we can write (compressible effects are 
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Figure 1 Time evolution of the diagnostic quantities for the fiducial set of runs ppm_hllc with different resolutions. The left panel 
shows the root mean square (RMS) Mach number, while the middle and left panels show, respectively, the ratios R xx /R yy and 
R xx /Rzz, R being the Reynolds stress tensor (equation 4). Since the z-direction is the anisotropic direction (it would play the role of 
the radial direction in a CCSN) the ratios R xx /R yy and R xx /R Z z, offer a global measure of the anisotropy of the flow at the largest 
scale. All of the quantities appear to have reached stationarity after time t > 3 and oscillate around their target values. All resolutions 
produce the same qualitative behavior. 
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Figure 4 Dissipation rate of the energy at the largest scales 
due to the turbulent cascade (not including direct dissipation by 
the numerical viscosity) as a function of resolution and for all of 
the schemes. The dissipation rate is normalized so as to be 1 
in the limit of large Reynolds numbers / resolution. At 128 3 
points all of the schemes show an error of less than 10%, with 
the HLLC schemes already close to the 2% level. 


negligible in our simulations): 

dtE(k) + d k II(k) ~ 0, (28) 

so that stationarity requires II(fc) ~ const. In particular, 
since energy is conserved, one finds II(A;) ~ (e). This 
means that, in the limit of large Reynolds numbers, the en¬ 
ergy decay rate depends only on the macroscopic properties 
of the flow (and in particular not on the nature of the vis¬ 
cosity), a fact that has also been verified numerically [59], 
The significance of this property and its importance for the 
modeling of turbulence cannot be overstated. 

In the context of CCSN simulations this means that the 
large scale kinetic energy, a crucial quantity for the dynam¬ 
ics of the explosion [15], can be faithfully captured even 
with simulations achieving modest Reynolds numbers. 

For an ILES, a basic requirement, then, is that a suffi¬ 
ciently high resolution should be achieved to correctly rep¬ 
resent the energy cascade at the largest scales. What qual¬ 
ifies as a sufficiently high resolution is of course depen¬ 
dent on the details of the closure built into the scheme (and 
on the accuracy required for the particular application). To 
quantify this, we can estimate the level of accuracy that can 
be reached at any given resolution, using our local simu¬ 
lations. In particular, we can study directly the energy flux 
across scales, defined by equation (16). This is shown in 
Figure 3 for all of the different runs,. 

As discussed before, we expect that II(fc) ~ (e) over an 
extended region in Fourier space should be a direct indica¬ 
tion that a simulation has been able to recover an inertial 
range. Perhaps not surprisingly, in light of previous results 
[24], we find that regions where II ~ (e) as wide as a few 
wave numbers 4 < k < 10 only appear at the highest reso¬ 
lutions (we will discuss the inertial range in more detail in 
Section 4.2). However, the amount of energy decaying from 
the largest scales reaches an asymptotic value much quicker 
than that implying that the total kinetic energy budget at the 
largest scales is well resolved even at modest resolutions. 


We can make a more quantitative statement concerning 
the energy decay rate by looking at the peak of the energy 
flux as a function of resolution, as shown in Figure 4. We 
can see that at 128 3 points all of the simulations have a de¬ 
viation from the asymptotic energy decay rate of less than 
10%. The least dissipative methods already have an error 
close to the 2% level. A comparison between PPM_HLLE 
and PPM_HLLC reveals the profound impact that the choice 
of the Riemann solver has even at relatively large scale 
(more on the dissipative properties of the different schemes 
in Section 4.3). 

4.2 Energy spectra 

Obviously, not all of the dynamics of turbulence can be 
reduced to the rate at which kinetic energy decays from the 
injection scale. The internal dynamics of the energy cas¬ 
cade, far from the injection scale and far from the dissipa¬ 
tion range, can also play an important role in many appli- 
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Figure 2 Square root of the magnitude of the vorticity, y/\V x v|, for four of the simulations with 512 3 resolution in a slice through the 
middle of the x-z plane at the final time of the simulations (t = 100). The panels show simulations using ppm_hlle_n 512 , 
ppm_hllc_n 512 , wenoz_hllc_n 5 12 , and tvd_hllc JJ 512 clockwise from the top left. The direction of the anisotropic driving is up in 
these figures. The colorcode goes linearly from 0 (no vorticity; dark colors) to 15 (light colors) and it is the same for all panels. 


cations. To analyze this aspect we consider in Figure 5 the 
energy spectrum of the velocity defined by equation ( 10 ). 
The spectra are compensated by A : 5 / 3 to highlight regions 
with Kolmogorov scaling, which might be expected in the 
inertial range. Since we want to focus on quantities that 
do not depend (or depend weakly) on the nature of the en¬ 
ergy injection at large scale, we show all of the spectra as 
a function of a dimensionless wave number, 512 k Ax. The 
rationale behind this normalization is that, first of all, we 
assume the Kolmogorov scale rj to be proportional to the 
grid spacing. Secondly, the 512 factor is introduced to have 
the dimensionless k, 512 /,- Ax coincide with the dimen¬ 
sional one for the highest resolution runs. With this choice, 
512 fc Ax = 512 corresponds to a wavelength of a single 


grid point, 512 fc Ax = 256 corresponds to a wavelength 
of two grid points and so on. 

Looking at any of the groups of runs in Figure 5, one 
can immediately notice that the spectra obtained at differ¬ 
ent resolutions do not collapse into a single curve in the 
dissipation region, as would be required by Kolmogorov’s 
first similarity hypothesis [50] ( cf, [60]). This lack of con¬ 
vergence in the dissipation region could be due to the non¬ 
linear viscosity of HRSC schemes. This, in turn, could re¬ 
sult in an anomalous scaling of 77 with the grid spacing. 
Such scaling has been reported in the past for ILES, but it 
is not very well understood [52]. The good agreement be¬ 
tween the three different groups of simulations employing 
the HLLC Riemann solver seems to support this hypothesis 
and suggests that the nonlinear viscosity introduced by the 
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Figure 3 Energy flux, as defined by equation (16), obtained with different numerical methods and resolutions. The energy flux is 
shown normalized to the average dissipation rate given by equation (24). From left to right and from top to bottom we show the results 
obtained with ppm.hllc, ppm_hllc_cflo . 8, ppm_hlle, tvd.hlle and wenoz_hllc. The bottom right panel show a comparison of 
all of the methods at 512 3 . All of the schemes show a good level of accuracy in the energy flux from the largest scales, with errors 
smaller than a few % already at low resolutions. The differences between the schemes become more marked at large wave numbers 
where the numerical dissipation starts to interfere with the energy cascade. 


Riemann solver is an important ingredient in setting this 
scaling. 

Convergence appears to be recovered at larger scales 
> 8 Ax (512fcAx < 64), but the spectra appear to be 
dominated by the bottleneck effect. This manifests itself 
as a bump in the compensated spectra extending from the 
dissipation range until the end of the inertial range, for 
the simulations that show one (e.g., until 512 A' Ax = 10 
for the HLLC runs), or until the energy injection scale 
(hi2 A: Ax = 4), for the simulations that show no or lit¬ 
tle inertial range (TVD.HLLE). The bottleneck effect is a 
viscous phenomenon which is also observed in direct nu¬ 
merical simulations. However, in the present context where 
viscosity is of numerical origin, it is at the very least ques¬ 
tionable if a pronounced bottleneck is a desirable feature of 
the modeling. In astrophysical flows, where the Reynolds 
numbers are typically very large, this pile up of energy at 
large scales is unphysical and could affect the quantitative 
and qualitative outcome of a simulation [13]. A quantifica¬ 
tion of the bottleneck effect in terms of the energy budget 
is discussed in Section 4.4. 

At even larger scales, an inertial range (E ~ k~ 5 / 3 and 
n ~ const, see Figure 3) seems to be recovered by the 
least dissipative schemes (PPM and WENOZ with HLLC) 
in the region 4 < k < 10. PPM.HLLE and TVD.HLLE 
have a more limited region, a few wave numbers at most, 
that could be interpreted as being an inertial range. We note 
that this resolution is not particularly high in comparison 
with state of the art DNS [59, 61], but it would already 


correspond to an extremely high resolution in global CCSN 
simulations that typically have ~ 30 — 70 zones across the 
turbulent region [13]. 

The overall behavior of the spectra, as obtained by all 
schemes, is consistent with Kolmogorov’s theory of turbu¬ 
lence. The anisotropic contributions to the angle-integrated 
spectra are too small to be detected in our data. 

4.3 Numerical viscosity 

At very small scales (~ several grid points) the dynam¬ 
ics is dominated by the numerical viscosity. This can be 
estimated from the residual of the energy equation (17) 
or, equivalently, by the effective numerical viscosity v{k) 
(equation 19). The latter is shown in Figure 6 for all 
schemes and resolutions. 

The first thing to notice is that the numerical viscosity 
provided by all numerical schemes is not constant, but dif¬ 
fers by roughly an order of magnitude between low and 
high k. Having a wave number dependent viscosity is a de¬ 
sirable feature expected in any LES model (explicit or oth¬ 
erwise). Nevertheless, this makes the definition and calcu¬ 
lation of the effective Reynolds number achieved in a simu¬ 
lation ambiguous. Meaningful ways to estimate it for ILES 
have been proposed [53] and they can be used to ease the 
comparison between different simulations and assess their 
quality. However, one has to be very careful while using 
any quoted “Reynolds number” from an ILES, to estimate 
things like the dynamical range achieved by a simulation, 
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Figure 5 Energy spectra (equation 10) obtained with different numerical methods and resolutions. The energy spectra are 
compensated by a k 5 / 3 spectrum, so that any region with Kolmogorov scaling should appear roughly flat. Furthermore, the spectra 
are all plotted as a function of the dimensionless wave number 512 k Ax (the 512 factor is introduced to have the dimensionless wave 
number coincide with the dimensional one for the 512 3 runs). The first five panels show the ppm.hllc (upper left), ppm_hllc_cflo . 8 
(upper center), ppm.hlle (upper right), tvd.hlle (lower left) and wenoz.hllc (lower center) group of runs. The last panel (lower 
right) shows a comparison of all of the methods at the highest resolution (512 3 ). An inertial range seems to be recovered only at the 
highest resolutions (perhaps with the exception of tvd_hlle where no inertial range is visible). All schemes employing the HLLC 
Riemann solver are in very good agreement. 


because the dissipative properties of ILES differ consider¬ 
ably from the ones of the true Navier-Stokes equations. 

Two other features can be observed in most of the numer¬ 
ical viscosity profiles. First, many of them exhibit a sudden 
reversal at high wave numbers. This is due to the fact that 
the numerical viscosity does not behave like a shear vis¬ 
cosity so that, although the numerical diffusion is strong at 
those scales, the numerical viscosity appears small because 
of a partial decoupling between vorticity and dissipation. 
Second, at high resolution and at the largest scales, the nu¬ 
merical viscosity is close to zero or even slightly negative. 
The reason is that the residual of equation (9) oscillates 
around zero and it is too small to be reliably extracted from 
our data: a much longer integration time would be needed 
to accumulate enough statistics for it. 

Finally, a comparison between the numerical viscos¬ 
ity reveals two interesting effects. First, by comparing 
PPM.HLLC and PPM.HLLE, we see that the choice of the 
Riemann solver affects the viscosity at basically all scales. 
Second, if we compare PPM.HLLC, PPM.HLLC.CFLO . 8 
and WENOZ.HLLC, we see that doubling the timestep ap¬ 
pears to have an effect comparable to the difference be¬ 
tween the PPM and WENOZ reconstructions at intermedi¬ 
ate scales (40 < k < 100). 

4.4 The energy distribution 

So far we have been concerned with the energy decay 
rate from the largest scales, which we have shown to be 


well captured by the IFES (Section 4.1), and with the en¬ 
ergy transfer in the inertial range, which we have seen to be 
described accurately only at much higher resolutions (Sec¬ 
tion 4.2). In a turbulent flow both of these aspects are im¬ 
portant and a good IFES should display a distribution of 
energy across vortical structures at different scales that is 
as close as possible to the asymptotic one. Obviously, there 
is a limit to the accuracy that any IFES can achieve at a 
fixed resolution. Here, we make this statement more quan¬ 
titative by considering the amount of kinetic energy that is 
well resolved by each simulation at a given resolution. 

We introduce the cumulative energy spectrum, the inte¬ 
gral of the energy spectrum: 

£(k)=[ k E( £)d£. (29) 

Jo 

This quantity is plotted in Figure 7, where it is normalized 
by 


+ OO 

E(k) d k (30) 

to obtain the cumulative distribution function of the kinetic 
energy. As a reference, we also show the cumulative energy 
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Figure 6 Numerical viscosity as a function of the wave number measured for all schemes and resolutions. The numerical viscosity is 
estimated using the procedure outlined in Section 2 and it is defined by equation (19). The different panels are, from left to right and 
from top to bottom, the results obtained with ppm_hllc, ppm_hllc_cflo . 8, ppm_hlle, tvd_hlle and wenoz_hllc. The bottom 
right panel show a comparison of all of the methods at 512 3 . The numerical viscosity shows large variations across the wave number 
space. The choice of the Riemann solver plays a role that is at least as important as the choice of the reconstruction method in 
affecting the numerical viscosity throughout the entire the spectrum. 


spectrum estimated from Kolmogorov’s theory: 

£k 4 ! (k)= [ £ K41 (£) d£, (31) 

Jo 

-E l K4l(fc) = 

-EpPM_HLLC-N512(fc) , if k < 4 , (32) 

£pPM_HLLC_N512(4) (j) , k > 4 . 

We find that as the resolution increases, all schemes ap¬ 
pear to be converging to the predictions of Kolmogorov’s 
theory. The results at finite resolution, however, are not en¬ 
couraging: at 64 3 only ~ 50% or less of the kinetic energy 
is in well resolved structures, while the other ~ 50% have 
piled up at rather large scale, with a cumulative excess of 
~ 20% at the grid scale, mostly because of the bottleneck 
effect. At higher resolutions, the amount of kinetic energy 
well captured by the ILES increases, but at 512 3 this is still 
only about 80% of the energy and there is still a cumulative 
excess of over ~ 5% at the grid scale {t ~ Ax). 

5 The 4/5-law 

The 4/5—law (equation 25) is not a-priori valid in the 
regime of turbulence we are considering. However, the 
4/5—law has been numerically verified to hold also in 
some situations outside the domain of validity of its deriva¬ 
tion. For instance, for isotropic mildly compressible de¬ 
caying [62] and driven [63] turbulence. In the anisotropic 


case, however, anisotropic contributions cannot be ex¬ 
cluded [64], although they are known to be subdominant 
in some important cases [65, 66, 67]. In this section we 
show that equation (25) is consistent with our data over a 
wide range of scales. 

We compute the 3 rd order structure functions of the ve¬ 
locity, defined by equation (23), in a rather simple way us¬ 
ing a random sample of 20,000 points in each of the 380 
3D data dumps of our simulations. At each time, we com¬ 
pute the 3 rd power of the velocity increments for each pair 
of points and accumulate and average in time the results 
in bins of size Ax’. The resulting structure functions are 
shown in Figure 8, compensated by — |r _1 (e) _1 , so that 
the resulting quantity should be equal to one if the 4/5—law 
is satisfied in our data. As was the case for the energy spec¬ 
tra, we assume r) ~ Ax and plot the structure functions 
versus r/ Ax. 

The degree with which the 4/5—law is satisfied in our 
data is very good. We see that anisotropic contributions 
only play a minor role in the angle-integrated formulation 
of the 4/5—law. This is in agreement with the incompress¬ 
ible DNS of [67] and has been known to be true also for 
Rayleigh-Benard convection in most regimes [68]. Our re¬ 
sults provide an important new example where this appears 
to hold true. Secondly, for all of our simulations at 512 3 , 
we find 

max|-^r _1 (e) _1 (<5u 3 )j =0 | (33) 
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Figure 7 Cumulated energy distribution (equation 29) for all methods and resolutions, normalized by a factor 2 /^ ms to be equal to 1 
for large k. As a reference for comparison we also plot the asymptotic profile expected from Kolmogorov’s theory (equation 31). The 
different panels are, from left to right and from top to bottom, the results obtained with ppnlhllc, ppm_hllc_cflo . 8, ppjylhlle, 
tvd.hlle and wenoz.hllc. The bottom right panel show a comparison of all of the methods at 512 3 . At low resolution all of the 
schemes show an excess of energy at intermediate scales, due to the bottleneck. Only at the highest resolution at least, roughly, 80% 
of the energy is correctly resolved. 


within 5% of 1. This level of accuracy is reached in 
DNS simulations achieving at least a Taylor micro-scale 
Reynolds number [67] 

u! A 

R x - - 300, (34) 

v 

where v is the kinematic viscosity, v! = ^u rms and 
A = (15 vu' 2 /(e)) 1 / 2 is the Taylor micro-scale. This corre¬ 
sponds to a large-scale Reynold numbers R = ~ R\, 

L = 1 being the domain size, in excess of —85,000. Reach¬ 
ing these Reynolds numbers in a DNS requires resolutions 
between 512 3 and 1024 3 using pseudo-spectral methods 
[69], This large-scale estimate of the Reynolds number is 
consistent with previous findings [53], although it is sev¬ 
eral orders of magnitude larger than the one that could be 
naively estimated using i/ max . For instance, for PPM_HLLC 
at 512 3 , zz max — 1.5 x 10 -3 and v rms — 0.4 giving 

v! L 

-- 150. (35) 

^max 

This apparent discrepancy is due to the fact that an ILES is 
not a DNS. As a consequence, different quantities that in 
a DNS depend on the Reynolds number, such as the dissi¬ 
pation rate or the Kolmogorov scale, behave as though the 
simulation had multiple values of the Reynolds number. 



512 k Aac 


Figure 9 3D (equation 10, blue) and 2D, transverse (equation 
26, red) energy spectra for the ppm_hllc simulations. The 
energy spectra are compensated by fc 5 / 3 to highlight eventual 
regions with Kolmogorov scaling. The spectra are plotted as a 
function of the dimensionless wave number 512 k Ax, as in 
Figure 5. Although E(k) and E ± (k ± ) have similar trends, the 
use of the transverse spectrum can overestimate the width of 
the bottleneck region. 


6 The transverse spectrum 

Finally, we want to comment on the use of 2D transverse 
spectra in 3D simulations, a practice typically employed in 
the analysis of turbulence in CCSN simulations [55, 23,47, 
13]. 

Figure 9 shows a comparison of the 2D transverse spec¬ 
trum E±(k±) from equation (26) and the 3D energy spec- 
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Figure 8 Compensated 3 rd order structure functions (equation 23) for all the numerical methods and resolutions. The structure 
functions are compensated and scaled so that they should be close to one where the 4/5-law (equation 25) is verified. The data is 
plotted as a function of the dimensionless separation r/ Ax. The first five panels show the ppnlhllc (upper left), ppm_hllc_cflo . 8 
(upper center), ppm.hlle (upper right), tvd_hlle (lower left) and wenoz.hllc (lower center) group of runs. The last panel (lower 
right) shows a comparison of all of the methods at the highest resolution (512 3 ). The 4/5-law is very well verified in our data 
suggesting that 1) anisotropic corrections are subdominant and 2) all of the simulations behave in a way consistent with large 
Reynolds numbers turbulence at the largest scales. 


trum from equation (10) for the PPNLHLLC simulations. 
The other runs (not shown) have the same qualitative be¬ 
havior. We can see that the transverse spectrum follows 
qualitatively the same trend as the 3D spectrum in terms 
of convergence. They are both roughly compatible with a 
Kolmogorov scaling, but the bottleneck appears to be more 
pronounced in the 2D spectrum than in the 3D spectrum. 
In particular, E±(k±) only shows a very small region that 
suggests an inertial range, 3 < k < 5 (as opposed to 
5 < k < 10 in E(k)). 

[13] concluded, also based on the analysis of 2D spec¬ 
tra, that turbulence in CCSN simulations is dominated by 
the bottleneck effect. Given the resolutions used in CCSN 
studies, our work supports their conclusion. However, in 
the light of Figure 9, we recommend that future studies sup¬ 
plement the analysis of 2D spectra with 3 rd order structure 
functions, that, as we have shown, can give a more accurate 
description of the energy cascade. 

7 Conclusions 

The details of the explosion mechanism of CCSNe have 
eluded our comprehension in spite of more than 50 years of 
studies [2, 3, 4, 5], Recent numerical advances [14, 48, 15, 
49] suggest that turbulence might play a fundamental role 
in tipping over the balance of the forces and lead to success¬ 
ful explosions (see also Appendix A). At the same time, the 
level of accuracy of current simulations, which employ the 
ILES methodology, is unclear [13], Turbulence in CCSNe 


is mildly compressible, but strongly anisotropic [14, 15], 
Simulations use rather dissipative numerical schemes (be¬ 
cause they have to deal with strong shock waves and com¬ 
plex microphysics) and relatively low resolution, a com¬ 
bination (anisotropic turbulence and dissipative schemes) 
that has not been systematically studied before. 

With the goal of assessing the reliability of ILES em¬ 
ployed in the study of CCSNe, as well as in other areas 
of physics and astrophysics, we performed a series of lo¬ 
cal simulations of driven, anisotropic, weakly compressible 
turbulence. We compared five commonly employed numer¬ 
ical schemes with different reconstruction methods, Rie- 
mann solvers, and time step size. Each was run at 4 dif¬ 
ferent resolutions ranging from 64 3 to 512 3 . Our analysis 
focused on the fidelity with which the turbulent cascade 
is represented in each model. In particular, we performed 
an analysis both in Fourier space (with the velocity power- 
spectra and the energy flux) and in physical space (with the 
3 rd order structure functions). Finally, we measured the nu¬ 
merical viscosity of each scheme from the residual of the 
specific kinetic energy equation. 

We found that, on the one hand, all of the numerical se¬ 
tups are able to accurately capture the decay rate of kinetic 
energy from the injection scale, with errors at the few % 
level already at 128 3 ( e.g ~ 2.5% for PPM_HLLC_N12 8). 
On the other hand, a large fraction of the energy is at unre¬ 
solved scales where it piles up due to the bottleneck effect 
and an inertial range appears only at the highest resolutions 
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(512 3 ). Even at this resolution, which would be difficult to 
achieve in global simulations, only roughly ~ 80% (the 
exact number depends on the scheme, see Section 4.4) of 
the energy is resolved, the remaining ~ 20% accumulates 
as excess energy at intermediate scales (the cumulative en¬ 
ergy excess at the grid scale alone is as large as ~ 5% of 
the total energy). 

Current CCSN simulations have resolutions of at most 
of 30 — 70 points covering the gain region [13] (the en¬ 
ergy injection scale). Based on our analysis we expect that 
at these resolutions even the energy decay rate from the 
largest scales will not be completely converged, but will 
show errors of up to tens of percent, depending on the nu¬ 
merical scheme (see Section 4.1). At smaller scales, the dy¬ 
namics is going to be completely dominated by the bottle¬ 
neck effect. This is in agreement with the findings of [13], 
based on the use of global simulations reaching a maximum 
resolution of 66 grid points covering radially the extent of 
the gain region. 

Based on our findings, we expect that, if the resolution in 
global simulations is increased by a factor ~ 2 from the one 
of [13], the decay rate will be converged to within a few % 
of the asymptotic value. This implies that the ratio between 
thermal and kinetic energy, a crucial quantity for the onset 
of the explosion, will also be converged to within a few %, 
at least when the energy injection rate changes slowly com¬ 
pared to the eddy turnover time (which is roughly ~ 20 ms 
in a CCSN [70, 23]). Unfortunately, while the lead up to 
explosion occurs over a larger timescale of a few hundred 
milliseconds, the transition to explosion can happen over 
much shorter timescales (one turnover time or less) [15], 
This means that the dynamics of the cascade over smaller 
time and length scales in the gain region also needs to be 
captured correctly since changes in the energy input rate on 
such short time scales will yield an inaccurate representa¬ 
tion of the energy on large scales due to the bottleneck ef¬ 
fect. This could require an increase of resolution by a factor 
~ 4 — 8 with respect to current high-resolution simulations. 
Additional work using semi-global or global simulations 
will be required to more firmly establish the resolutions re¬ 
quirements at the transition of the explosion. 

Concerning the properties of anisotropic turbulence in 
our simulations, we found anisotropy contributions to the 
energy spectrum and to the angle-averaged formulation of 
the 4/5—law to be subdominant: the accuracy with which 
the 4/5—law is satisfied is limited only by the employed 
resolution and the energy spectrum appears to be consis¬ 
tent with Kolmogorov k~ 5 ' 3 scaling. We also found the 
transverse energy spectrum with respect to the direction of 
anisotropy, a quantity typically computed in CCSN simu¬ 
lations, to overestimate the bottleneck with respect to the 
angle-integrated 3D spectrum. For this reason, we recom¬ 
mend future studies of CCSN to supplement (or replace) 
the analysis of the transverse spectrum with the analysis 


of the 3 rd order, angle-integrated, structure function (or, 
where possible, with the 3D spectrum itself). 

Our results are, of course, dependent on the choice of the 
numerical scheme. In particular, we found significant dif¬ 
ferences in the dissipative properties of schemes employing 
the HLLC Riemann solver with respect to schemes using 
the more dissipative HLLE solver. The reconstruction order 
of the scheme is also important, although, while significant 
differences are found between TVD and PPM, the differ¬ 
ences between PPM and WENOZ are much more minute 
(despite WENOZ being significantly more computationally 
expensive than PPM). In the end, none of the schemes we 
considered seems to be able to yield an accurate representa¬ 
tion of the kinetic energy distribution across different scales 
at an affordable resolution for global CCSN simulations. 
A possible way forward would be to adopt low-dissipation 
numerical schemes especially designed for the use in ILES, 
such as the methods proposed by [71, 72] or [73], Imple¬ 
menting and testing these schemes will be subject of future 
work. 

An important limitation of the present work is that we 
considered a very idealized setup. On the one hand, this al¬ 
lows us to benchmark the behavior of ILES in a controlled 
environment. On the other hand, our simulations cannot 
fully capture all features of the turbulent convective flow 
in a CCSN. Unlike the situation in a CCSN, our local sim¬ 
ulations did not include a vertical advective velocity field 
that is due to the accretion of the stellar mantle. However, 
the advective velocities are nearly constant in the regions 
of interest and Galilean invariance ensures that our results 
are unaffected. More limiting is the local nature of our sim¬ 
ulations and the inevitable choice of boundary conditions. 
Moreover, our simulations could not take into account spa¬ 
tial variations in gravity and the large-scale radial conver¬ 
gence of the flow in globally spherical problems like col¬ 
lapsing stars. Addressing these issues will also be subject 
of future work. 

Appendix A: The role of turbulence in 
core-collapse supernova 
explosions 

In this appendix we present a brief discussion of the im¬ 
portance of turbulent pressure in the explosion of massive 
stars. To set the stage, we will briefly summarize what 
is known of the dynamics of the most common class of 
CCSNe that are relevant for our later discussion. This is 
done for the benefit of readers that are not supernova spe¬ 
cialists and it is not meant to be a complete review of the 
status of the field, for which we refer, instead, to the reviews 
of [3] and [4], Next, we will discuss the role of turbulence 
and, in particular, of turbulent pressure on the explosion 
mechanism, in light of some recent results [14, 15], 
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A.1 The neutrino mechanism 

Towards the end of their evolution, massive stars form mas¬ 
sive (~ 1.5 solar mass) iron cores at their center. Since the 
iron nucleus has the largest binding energy per nucleon, 
no energy can be extracted from nuclear fusion beyond 
iron. The iron core is essentially inert and supported against 
gravity only by the degeneracy pressure of relativistic elec¬ 
trons. The mass of the iron core increases with time as 
more iron-group material is added by silicon shell burning. 
Electron capture on protons, which becomes energetically 
favorable at high densities, depletes the core of electrons 
and thus reduces the pressure supporting it against grav¬ 
ity. Eventually, the core becomes dynamically unstable and 
collapses. 

During the collapse, the subsonically collapsing inner 
core (~ 0.5 solar masses) contracts until it reaches den¬ 
sities comparable to that in atomic nuclei (~ 4 — 5 x 
10 14 g cm -3 ). At this point, the nuclear equation of state 
stiffens (due to the strong nuclear force). This halts the col¬ 
lapse of the inner core. It stops, bounces back and a proto¬ 
neutron star (PNS) is formed. The outer core, however, is 
still collapsing supersonically and a strong shock wave is 
launched at the interface between the inner and outer core. 

It was once thought that this shock wave would travel 
outwards dynamically, depositing its energy in the outer 
layers of the star, causing the explosion. However, multi¬ 
ple numerical simulations performed over multiple decades 
have consistently shown that the initial shock fails to ex¬ 
plode the star. Instead, it stalls due to energy losses to the 
dissociation of heavy nuclei into free nucleons and to the 
emission of neutrinos that stream away from the neutrino- 
semitransparent regions behind the shock [74], The shock 
generally stalls within only a few tens of milliseconds of 
core bounce and turns into an accretion shock standing at a 
radius of ~ 100 — 200 km. The accretion rate through the 
shock is so high (a fraction of a solar mass per second) that, 
if nothing revitalizes the shock within ~ 1 — 2 seconds, the 
gravitational force would overwhelm the nuclear repulsion 
force, collapsing the core of the supernova to a black hole, 
precluding explosion (e.g., [75]). 

During this time, however, the PNS will release a signif¬ 
icant fraction of its binding energy in the form of neutrinos 
(of order 10 46 J). Converting a few percent of that energy 
into kinetic energy would be enough to unbind the stellar 
envelope and power the supernova explosion. In the stan¬ 
dard neutrino mechanism it is theorized that a small frac¬ 
tion neutrinos emitted from the edge of the protoneutron 
star is re-absorbed in the region right behind the stalled 
shock. The deposition of neutrino energy leads to higher 
thermal pressure so that the shock can eventually overcome 
the ram pressure of accretion and accelerates in a run-away 
process [74, 76]. Turbulence in the heating region behind 
the shock increases the time a fluid parcel spends in that 
region and, importantly, turbulent pressure helps in over¬ 
coming the ram pressure of accretion (see next Section and 


[23]). It is, however, presently unclear if neutrino heating 
(even if aided by turbulence in launching the explosion) is 
able to provide enough energy to power the explosions to 
the energies inferred from astronomical observations. 

A.2 Turbulent pressure and the Rankine-Hugoniot 
conditions 

Simulations [8, 14, 15] have shown that turbulence and, in 
particular, turbulent pressure behind the shock, could play 
an important role in aiding the explosion. To see why this 
is the case, we consider the Rankine-Hugoniot momentum 
condition for a standing accretion shock in a supernova 
core, 

s[PdV r d - PuV r u ] = Pd(y r d ) 2 +p d - pu(v r u ) 2 - p u , (36) 

where s is the shock speed and - d and ■„ denote the down¬ 
stream and upstream values respectively. For the purpose 
of our discussion, we can assume the upstream gas to be 
cold and free-falling: 

, GM 

Pu = o K) 2 = — • (37) 

For the shock to expand we must then have 

PdKf^Pu™. (38) 

r 

In the presence of turbulence, [14] suggested to modify 
equation (38) in a way akin to a Reynolds decomposition 
and write it as 

Pd{v r df + p d {Sv r d ) 2 +p d > pj^~, (39) 

r 

where v is the average velocity and <5v = v—v is the turbu¬ 
lent velocity. Although not entirely rigorous, equation (39) 
has been shown to be well verified in the numerical simu¬ 
lations if angular averages are used to compute the respec¬ 
tive quantities [14, 15]. [15] have shown that the turbulent 
pressure expressed in this fashion can exceed 50% of the 
thermal pressure, making a very significant contribution to 
the momentum balance in (39). 

Going beyond the arguments of [14], we can reinterpret 
equation (39) as being the Rankine-Hugoniot condition for 
a fluid with a modified equation of state, which has two 
separate internal degrees of freedom: thermodynamical and 
turbulent. To this aim, we express 5v d in terms of the spe¬ 
cific turbulent energy 

eturb = ^l<5v| 2 (40) 

noting that 

|<5v| 2 := (5v r ) 2 + (Sv 9 ) 2 + (6v*) 2 , 


(41) 
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and using the fact that 

0 5v r ) 2 ~ (Sv 6 ) 2 + (Sv 4 ’) 2 (42) 

in CCSN turbulence, to obtain 

(5v r ) 2 ~ ^|<5v| 2 = e t urb- (43) 

Assuming the pressure varies like p = (j — l)pe, and sub¬ 
stituting (43) into (39), we find 

PdiVd) + (7th ^)pd&d 4“ ('Yturb ^)Pd &turb ^ Pu 5 

r 

(44) 

where 7 th — 4/3 is the thermodynamical adiabatic index, 
e,i is the downstream thermal energy and 7t U rb = 2 is 
the equivalent adiabatic index associated with anisotropic 
CCSN turbulence. Since 7 t ur b > 7 th, we see that turbulent 
energy is more efficient, per unit specific internal energy, at 
pushing the shock than thermal energy. 

We point out that, if equation (42) is dropped and tur¬ 
bulence is assumed to be isotropic, then 7t m b = 5/3, 
which is still larger than 7 t h, but not as large as for the 
anisotropic case. This is a simple consequence of the fact 
that anisotropic turbulence has an anisotropic pressure, 
which is stronger in the radial direction. 

In both cases, since the total energy is conserved, the rel¬ 
evant quantity is the ratio eturb /e. From standard turbulent 
phenomenology we expect that this ratio will only depend 
on macroscopic parameters, such as the net heating rate, the 
accretion rate and so on, and not on the details of the vis¬ 
cosity. For this reason, we expect this ratio to be correctly 
captured in ILES achieving a sufficiently high resolution. 

As a final remark, we point out that a similar argument 
has been recently proposed by [49] who formulated their 
equations in terms of the turbulent Mach number, as op¬ 
posed to the turbulent energy. 

Competing interests 

The authors declare that they have no competing interests. 

Author’s contributions 

DR ran the simulations, performed the analysis of the data and wrote the 
basic draft of this paper. SMC assisted with the use of the FLASH code. 
CDO had the original idea that started this investigation. All of the authors 
took part in discussions concerning the results and contributed corrections 
and improvements on the early draft of the manuscript. 

Acknowledgements 

We acknowledge helpful discussions with E. Abdikamalov, W. D. Arnett, 

A. Burrows, R. Fisher, C. Meakin, R Mosta, J. Murphy, M. Norman, and 
L. Roberts. This research was partially supported by the National Science 
Foundation under award nos. AST-1212170 and PHY-1151197 and by the 
Sherman Fairchild Foundation. The simulations were performed on the 
Caltech compute cluster Zwicky (NSF MRI-R2 award no. PHY-0960291), on 
the NSF XSEDE network under allocation TG-PHY100033, and on 
NSF/NCSA BlueWaters under NSF PRAC award no. ACI-1440083. The 
software used in this work was in part developed by the DOE NNSA-ASC 
OASCR Flash Center at the University of Chicago. 


References 

1. Clausen, D., Piro, A.L., Ott, C.D.: The Black Hole Formation Probability. 
ApJ 799, 190 (2015). doi:10.1088/0004-637X/799/2/190. 1406.4869 

2. Woosley, S.E., Janka, H.-T.: The physics of core-collapse supernovae. 
Nature Physics 1, 147 (2005). astro-ph/0601261 

3. Janka, H.-T., Hanke, F., Htidepohl, L., Marek, A., Muller, B., 
Obergaulinger, M.: Core-collapse supernovae: Reflections and 
directions. Prog. Th. Exp. Phys. 2012 ( 1 ), 010000 (2012). 1211.1378 

4. Burrows, A.: Colloquium: Perspectives on core-collapse supernova 
theory. Rev. Mod. Phys. 85 , 245 (2013). 

doi:1 0.1103/RevModPhys.85.245. 1210.4921 

5. Foglizzo, T., Kazeroni, R., Guilet, J., Masset, F., Gonzalez, M., Others: 
The explosion mechanism of core-collapse supernovae: progress in 
supernova theory and experiments (2015). 1501.01334 

6. Janka, H.-T.: Explosion Mechanisms of Core-Collapse Supernovae. 
Ann. Rev. Nuc. Par. Sci. 62 , 407 (2012). 

doi:1 0.1146/annurev-nucl-102711 -094901. 1206.2503 

7. Herant, M.: The convective engine paradigm for the supernova 
explosion mechanism and its consequences. Phys. Rep. 256 , 1 17 
(1995). doi:1 0.1016/0370-1573(94)00105-C 

8. Burrows, A., Hayes, J., Fryxell, B.A.: On the Nature of Core-Collapse 
Supernova Explosions. ApJ 450 , 830 (1995) 

9. Janka, H.-T., Muller, E.: Neutrino heating, convection, and the 
mechanism of Type-ll supernova explosions. A&A 306, 167 (1996) 

10. Foglizzo, T., Scheck, L., Janka, H.-T.: Neutrino-driven Convection 
versus Advection in Core-Collapse Supernovae. ApJ 652 , 1436 (2006). 
doi:10. 1086/508443 

11. Blondin, J.M., Mezzacappa, A., DeMarino, C.: Stability of Standing 
Accretion Shocks, with an Eye toward Core-Collapse Supernovae. ApJ 
584 , 971 (2003). doi: 1 0.1086/345812 

12. Foglizzo, T., Galletti, R, Scheck, L., Janka, H.-T.: Instability of a Stalled 
Accretion Shock: Evidence for the Advective-Acoustic Cycle. ApJ 654 , 
1006 (2007). doi: 10.1086/509612 

13. Abdikamalov, E., Ott, C.D., Radice, D., Roberts, L.F., Haas, R., 
Reisswig, C., Mosta, R, Klion, H., Schnetter, E.: Neutrino-driven 
Turbulent Convection and Standing Accretion Shock Instability in 
Three-dimensional Core-collapse Supernovae. ApJ 808, 70 (2015). 
doi:1 0.1088/0004-637X/808/1 /70. 1409.7078 

14. Murphy, J.W., Dolence, J.C., Burrows, A.: The Dominance of 
Neutrino-driven Convection in Core-collapse Supernovae. ApJ 771,52 
(2013). doi:1 0.1088/0004-637X/771/1/52. 1205.3491 

15. Couch, S.M., Ott, C.D.: The Role of Turbulence in Neutrino-driven 
Core-collapse Supernova Explosions. ApJ 799, 5 (2015). 

doi:1 0.1088/0004-637X/799/1/5. 1408.1399 

16. Gamier, E., Adams, N., Sagaut, R: Large Eddy Simulation for 
Compressible Flows. Springer (2000) 

17. Grinstein, F.F., Margolin, L.G., Rider, W.J.: Implicit Large Eddy 
Simulation: Computing Turbulent Fluid Dynamics, p. 578. Cambridge 
University Press (2011) 

18. Yakhot, V., Zakharov, V.: Hidden conservation laws in hydrodynamics; 
energy and dissipation rate fluctuation spectra in strong turbulence. 
Phys. D 64(4), 379(1993) 

19. She, Z., Jackson, E.: On the universal form of energy spectra in fully 
developed turbulence. Phys. Fluids 5(7), 1526 (1993) 

20. Falkovich, G.: Bottleneck phenomenon in developed turbulence. Phys. 
Fluids 6(4), 1411-1414(1994) 

21. Verma, M.K., Donzis, D.: Energy transfer and bottleneck effect in 
turbulence. J. Phys. A 40(16), 4401 (2007) 

22. Frisch, U., Kurien, S., Pandit, R., Pauls, W., Ray, S., Wirth, A., Zhu, 

J. -Z.: Hyperviscosity, Galerkin Truncation, and Bottlenecks in 
Turbulence. Phys. Rev. Lett. 101(14), 144501 (2008). 

doi:1 0.1103/PhysRevLett.l 01.144501 

23. Couch, S.M., O’Connor, E.P.: High-resolution Three-dimensional 
Simulations of Core-collapse Supernovae in Multiple Progenitors. ApJ 
785, 123 (2014). doi:10.1088/0004-637X/785/2/123. 1310.5728 

24. Sytine, I.V., Porter, D.H., Woodward, P.R., Hodson, S.W., Winkler, 

K. -H.: Convergence Tests for the Piecewise Parabolic Method and 
Navier-Stokes Solutions for Homogeneous Compressible Turbulence. 

J. Comp. Phys. 158, 225 (2000). doi :10.1006/jcph. 1999.6416 

25. Gamier, E., Mossi, M., Sagaut, P.: On the use of shock-capturing 
schemes for large-eddy simulation. Journal of Computational ... 311, 



Rad ice etal. 


Page 16 of 17 


273-311 (1999) 

26. Johnsen, E., Larsson, J., Bhagatwala, A.V., Cabot, W.H., Moin, P., 
Olson, B.J., Rawat, P.S., Shankar, S.K., Sjogreen, B., Yee, H.C., Zhong, 
X., Lele, S.K.: Assessment of high-resolution methods for numerical 
simulations of compressible turbulence with shock waves. Journal of 
Computational Physics 229 ( 4 ), 1213-1237 (2010). 

doi:1 0.1016/j.jcp.2009.10.028 

27. Arnett, D., Meakin, C., Young, P.A.: Turbulent Convection in Stellar 
Interiors. II. The Velocity Field. ApJ 690 , 1715 (2009) 

28. Casciola, C.M., Gualtieri, P., Jacob, B., Piva, R.: The residual 
anisotropy at small scales in high shear turbulence. Physics of Fluids 
19(10), 101704 (2007). doi:1 0.1063/1.2800043 

29. Colella, P., Woodward, P.R.: The Piecewise Parabolic Method (PPM) for 
Gas-Dynamical Simulations. J. Comp. Phys. 54, 174 (1984) 

30. Toro, E.F.: Riemann Solvers and Numerical Methods for Fluid 
Dynamics. Springer, Berlin (1999) 

31. Schmidt, W., Hillebrandt, W., Niemeyer, J.C.: Numerical dissipation and 
the bottleneck effect in simulations of compressible isotropic 
turbulence. Computers & Fluids 35(4), 353 (2006) 

32. Thornber, B., Mosedale, A., Drikakis, D.: On the implicit large eddy 
simulations of homogeneous decaying turbulence. Journal of 
Computational Physics 226 ( 2 ), 1902-1929 (2007). 

doi :1 0.1016/j.jcp.2007.06.030 

33. Fryxell, B., Olson, K., Ricker, R, Timmes, F.X., Zingale, M., Lamb, D.Q., 
MacNeice, R, Rosner, R., Truran, J.W., Tufo, H.: FLASH: An Adaptive 
Mesh Hydrodynamics Code for Modeling Astrophysical Thermonuclear 
Flashes. ApJS 131 , 273 (2000). doi: 1 0.1086/317361 

34. Dubey, A., Antypas, K., Ganapathy, M.K., Reid, L.B., Riley, K., Sheeler, 
D., Siegel, A., Weide, K.: Extensible component-based architecture for 
flash, a massively parallel, multiphysics simulation code. Parallel 
Comput. 35 , 512 (2009). doi:10.1016/j.parco.2009.08.001 

35. Lee, D., Tzeferacos, R, Couch, S.M., Bachan, J., Daley, C., Fatenejad, 
M., Flocke, N., Lamb, D., Weide, K., Dubey, A.: Flash: A multi-physics 
code for adaptive mesh computational fluid dynamics in astrophysics. 
To be submitted to ApJS(2014) 

36. Colella, P.: Multidimensional upwind methods for hyperbolic 
conservation laws. J. Comput. Phys. 87, 171-200 (1990). 
doi :1 0.1016/0021 -9991 (90)90233-Q 

37. Lee, D., Deane, A.E.: An unsplit staggered mesh scheme for 
multidimensional magnetohydrodynamics. J. Comput. Phys. 228, 
952-975 (2009). doi:1 0.1016/j.jcp.2008.08.026 

38. Lee, D.: A solution accurate, efficient and stable unsplit staggered 
mesh scheme for three dimensional magnetohydrodynamics. J. Comp. 
Phys. 243, 269 (2013). doi:1 0.1016/j.jcp.2013.02.049. 1303.6988 

39. Borges, R., Carmona, M., Costa, B., Don, W.S.: An improved weighted 
essentially non-oscillatory scheme for hyperbolic conservation laws. 
Journal of Computational Physics 227 , 3191-3211 (2008). 

doi:1 0.1016/j.jcp.2007.11.038 

40. Einfeldt, B.: On Godunov-Type Methods for Gas Dynamics. SIAM 
Journal on Numerical Analysis 25 , 294-318 (1988). 

doi:10.1 137/0725021 

41. Toro, E.F., Spruce, M., Speares, W.: Restoration of the contact surface 
in the HLL-Riemann solver. Shock Waves 4, 25-34 (1994). 

doi :1 0.1007/BF01414629 

42. Uhlenbeck, G., Ornstein, L.: On the theory of the brownian motion. 
Phys. Rev. 36 , 823-841 (1930). doi:1 0.1103/PhysRev.36.823 

43. Eswaran, V., Pope, S.: An examination of forcing in direct numerical 
simulations of turbulence. Computers & Fluids 16(3), 257-278 (1988) 

44. Federrath, C., Roman-Duval, J., Klessen, R.S., Schmidt, W., Mac Low, 
M.-M.: Comparing the statistics of interstellar turbulence in simulations 
and observations. Astronomy and Astrophysics 512 , 81 (2010). 

doi :1 0.1051 /0004-6361 /200912437 

45. Verhoeven, J., Wiesehofer, T., Stellmach, S.: Anelastic versus Fully 
Compressible Turbulent Rayleigh-Benard Convection. ApJ 805 , 62 
(2015). doi: 10.1088/0004-637X/805/1 /62. 1501.01237 

46. Murphy, J.W., Meakin, C.: A Global Turbulence Model for 
Neutrino-driven Convection in Core-collapse Supernovae. ApJ 742, 74 
(2011). doi: 10.1088/0004-637X/742/2/74. 1106.5496 

47. Handy, T., Plewa, T., Odrzywotek, A.: Toward Connecting 
Core-collapse Supernova Theory with Observations. I. Shock Revival 
in a 15 M 0 Blue Supergiant Progenitor with SN 1987A Energetics. ApJ 


783, 125 (2014). doi:10.1088/0004-637X/783/2/125 

48. Couch, S.M., Ott, C.D.: Revival of The Stalled Core-Collapse 
Supernova Shock Triggered by Precollapse Asphericity in the 
Progenitor Star. ApJ 778 , 7 (2013) 

49. Muller, B., Janka, H.-T.: Non-radial instabilities and progenitor 
asphericities in core-collapse supernovae. Mon. Not. R. Astron. Soc. 
448 , 2141-2174 (2015). doi:10.1093/mnras/stv101. 1409.4783 

50. Frisch, U.: Turbulence: The Legacy of A. N. Kolmogorov. Cambridge 
University Press (1996) 

51. Fureby, C., Grinstein, F.F.: Monotonically Integrated Large Eddy 
Simulation of Free Shear Flows. AIAA Journal 37(5), 544-556 (1999). 
doi:1 0.2514/2.772 

52. Aspden, A., Nikiforakis, N., Dalziel, S., Bell, J.: Analysis of implicit LES 
methods. CAMCS 3(1), 103 (2009) 

53. Zhou, Y., Grinstein, F.F., Wachtor, A.J., Haines, B.M.: Estimating the 
effective Reynolds number in implicit large-eddy simulation. Physical 
Review E 89(1), 013303 (2014). doi:10.1 103/PhysRevE.89.013303 

54. Biferale, L., Procaccia, I.: Anisotropy in turbulent flows and in turbulent 
transport. Physics Reports 414(2-3), 43-164 (2005). 
doi:10.1016/j.physrep.2005.04.001 

55. Dolence, J.C., Burrows, A., Murphy, J.W., Nordhaus, J.: Dimensional 
Dependence of the Hydrodynamics of Core-collapse Supernovae. ApJ 
765 , 110 (2013). 1210.5241 

56. Porter, D.H., Woodward, P.R.: High-resolution simulations of 
compressible convection using the piecewise-parabolic method. ApJS 
93, 309 (1994). doi :10.1086/192057 

57. Vincent, A., Meneguzzi, M.: The satial structure and statistical 
properties of homogeneous turbulence. Journal of Fluid Mechanics 
225 , 1 (1991). doi:1 0.1017/S0022112091001957 

58. Ishihara, T., Kaneda, Y., Yokokawa, M., Itakura, K., Uno, A.: Small-scale 
statistics in high-resolution direct numerical simulation of turbulence: 
Reynolds number dependence of one-point velocity gradient statistics. 
Journal of Fluid Mechanics 592, 335-366 (2007). 

doi:1 0.1017/S0022112007008531 

59. Kaneda, Y., Ishihara, T., Yokokawa, M., Itakura, K., Uno, A.: Energy 
dissipation rate and energy spectrum in high resolution direct 
numerical simulations of turbulence in a periodic box. Phys. Fluids 15 , 
21 (2003). doi:10.1063/1. 1539855 

60. Gotoh, T., Fukayama, D., Nakano, T.: Velocity field statistics in 
homogeneous steady turbulence obtained using a high-resolution 
direct numerical simulation. Physics of Fluids 14(3), 1065 (2002) 

61. Federrath, C.: On the universality of supersonic turbulence. Mon. Not. 

R. Astron. Soc. 436(2), 1245-1257 (2013). doi:10.1093/mnras/stt1644 

62. Porter, D., Pouquet, A., Woodward, R: Measures of intermittency in 
driven supersonic flows. Physical Review E 66(2), 1-12 (2002). 
doi:1 0.1103/PhysRevE.66.026301 

63. Benzi, R., Biferale, L.. Fisher, R.T., Kadanoff, L.P., Lamb, D.Q., Toschi, 
F.: Intermittency and Universality in Fully Developed Inviscid and 
Weakly Compressible Turbulent Flows. Phys. Rev. Lett. 100(23), 
234503 (2008). doi:1 0.1103/PhysRevLett.l00.234503. 0709.3073 

64. Biferale, L., Daumont, I., Lanotte, a., Toschi, F.: Anomalous and 
dimensional scaling in anisotropic turbulence. Physical Review E 66 ( 5 ), 
056306 (2002). doi:1 0.1103/PhysRevE.66.056306 

65. Calzavarini, E., Toschi, F., Tripiccione, R.: Evidences of 
bolgiano-obhukhov scaling in three-dimensional rayleigh-benard 
convection. Phys. Rev. E 66, 016304 (2002). 

doi:10.1 103/PhysRevE.66.016304 

66. Biferale, L., Calzavarini, E., Toschi, F., Tripiccione, R.: Universality of 
anisotropic fluctuations from numerical simulations of turbulent flows. 
Europhysics Letters (EPL) 64 ( 4 ), 461^t67 (2003). 

doi:1 0.1209/epl/i2003-00233-9 

67. Kaneda, Y., Yoshino, J., Ishihara, T.: Examination of Kolmogorov’s 4/5 
Law by High-Resolution Direct Numerical Simulation Data of 
Turbulence. Journal of the Physical Society of Japan 77(6), 064401 
(2008). doi:1 0.1143/JPSJ.77.064401 

68. Lohse, D., Xia, K.-Q.: Small-Scale Properties of Turbulent 
Rayleigh-Benard Convection. Annual Review of Fluid Mechanics 42(1), 
335-364 (2010). doi:10.1 146/annurev.fluid.010908.165152 

69. Donzis, D.a., Yeung, P.K., Sreenivasan, K.R.: Dissipation and enstrophy 
in isotropic turbulence: Resolution effects and scaling in direct 
numerical simulations. Phys. Fluids 20(4), 045108 (2008). 


Rad ice etal. 


Page 17 of 17 


doi:1 0.1063/1.2907227 

70. Ott, C.D., Abdikamalov, E., Mosta, P., Haas, R., Drasco, S., O'Connor, 
E.P., Reisswig, C., Meakin, C.A., Schnetter, E.: General-relativistic 
Simulations of Three-dimensional Core-collapse Supernovae. ApJ 768 , 
115 (2013). 1210.6674 

71 . Hickel, S., Adams, N.a., Domaradzki, J.A.: An adaptive local 
deconvolution method for implicit LES. Journal of Computational 
Physics 213 ( 1 ), 413-436 (2006). doi:1 0.1016/j.jcp.2005.08.01 7 

72. Martin, M.P., Taylor, E.M., Wu, M., Weirs, V.G.: A bandwidth-optimized 
WENO scheme for the effective direct numerical simulation of 
compressible turbulence. Journal of Computational Physics 220 ( 1 ), 
270-289 (2006). doi:1 0.1016/j.jcp.2006.05.009 

73. Thornber, B., Mosedale, a., Drikakis, D., Youngs, D., Williams, R.J.R.: 
An improved reconstruction method for compressible flows with low 
Mach number features. Journal of Computational Physics 227 ( 10 ), 
4873-4894 (2008). doi:1 0.1016/j.jcp.2008.01.036 

74. Bethe, H.A.: Supernova mechanisms. Rev. Mod. Phys. 62, 801 (1990) 

75. O’Connor, E., Ott, C.D.: Black Hole Formation in Failing Core-Collapse 
Supernovae. ApJ 730 , 70 (2011) 

76. Pejcha, O., Thompson, T.A.: The Physics of the Neutrino Mechanism of 
Core-collapse Supernovae. ApJ 746 , 106 (2012). 

doi :1 0.1088/0004-637X/746/1 /1 06. 1103.4864 


